
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE)    #  erfc=1-erf

rmse = function(sim, obs){
  sqrt(mean((sim - obs)^2))
}

sech= function(x) 1/(cosh(x))

BC=function(BC_sub,suction_p) {
  theta_r=BC_sub$theta_r
  theta_s=BC_sub$theta_s
  phi_s=BC_sub$phi_s
  n=BC_sub$n
  
  theta=NULL
  for (i in seq(length(suction_p)) ) {
    
    if (suction_p[i]>phi_s ) {  
      theta[i]=theta_r+(theta_s-theta_r)* (phi_s/abs(suction_p[i]))^n
    } else {
      theta[i] = theta_s
    }
  }
  
  return(theta)
}

BC2=function(BC_sub,suction_p) {
  theta_r=BC_sub$theta_r
  theta_s=BC_sub$theta_s
  psi_s=BC_sub$psi_s
  n=BC_sub$n
  
  theta=NULL
  for (i in seq(length(suction_p)) ) {
    
    if (suction_p[i]>psi_s ) {  
      theta[i]=theta_r+(theta_s-theta_r)* (abs(psi_s/suction_p[i]))^n
    } else {
      theta[i] = theta_s
    }
  }
  return(theta)
}

vG=function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = thr + (ths-thr)/(1+(alp*h)^n)^(1-1/n)
  return(t)
}

kom <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = thr + (ths-thr)*0.5* erfc((log(h)-alp)/(sqrt(2)*n))
  return(t)
}

lom <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = thr + (ths-thr)*0.5*(1-tanh((log(h)-alp)/n))
  return(t)
}

atm <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = thr + (ths-thr)/2 - (ths-thr)/pi * atan((log(h)-alp)/n)
  return(t)
}

gdm <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  gd  <- function(t) atan(sinh(t))
  t = thr + (ths-thr)/2 - (ths-thr)/pi *gd((log(h)-alp)/n)
  return(t)
}

# change constant coefficient pjli 2022/03/18
e2m <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = ths -(ths-thr) * (1 + exp(-(log(h)-alp)/n))^(-1)
  return(t)
}

sqm <- function(para,h) {
  thr = para$thr
  ths = para$ths
  alp = para$alp
  n   = para$n
  t = thr + (ths-thr)/2 - (ths-thr)/2 * (log(h)-alp)/n * (1 + ((log(h)-alp)/n)^2)^(-1/2)
  return(t)
}

lr <- function(para,sand,clay) {
  a1 = para$a1
  a2 = para$a2
  b = para$b
  t = sand*a1 + clay*a2 + b
  return(t)
}











